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Greenbelt, Maryland 


INTRODUCTION 


This document contains a complete description of the Brouwer-Lyddane Orbit 
generator routine (BROWRO) which was written as a part of the Definitive Orbit 
Determination System (DODS). The routine accepts as input a set of Brouwer 
Mean elements at a specified epoch and outputs the position and velocity vectors 
at requested times along with intermediate data if requested. The formulation 
for the routine was derived from the following sources: 

1. Brouwer, Dirk. "Solution of the Problem of Artificial Satellite Theory Without 
Drag." The Astronomical Journal , Vol. 64, (Oct. 1959), 378-397. 

2. Lyddane, R. H. "Small Eccentricities or Inclinations in the Brouwer Theory 

of Artificial Satellite." The Astronomical Journal, Vol. 68, (June 1963), 555- 
558. " 

3. Siry, J. W. "A Mathematical Formulation of the Brouwer-Lyddane Orbit 
Theory." 

In addition to the basic Brouwer Theory, the BROWRO routine calls sub- 
routines for computing Delta L Drag and complementary perturbations in order 
to allow for the inclusion of these effects on the satellite's orbit. This document 
contains a description of these subroutines. 

The Brouwer-Lyddane Orbit generator routine was checked out first in a 
stand-alone environment and then within the DODS system. DODS is a Fortran 
language system of programs that execute under MVT on the S/360 75 and 95 
computers. The size of the system made it necessary to structure it in segments 
that could overlay each other in the computer as functional needs changed during 
execution. Thus, instead of one complete program or set of programs residing 
in core storage during the entire running time, only those portions needed are in 
core. Because interfacing with DODS was a major consideration in the program- 
ming, some techniques are utilized that would not be necessary in a stand-alone 
environment. For example, in BROWRO many variables which could be assumed 
to contain valid data once they had been computed in a stand-alone environment 
had to be recomputed when running under DODS since there was a possibility of 
their being destroyed by the DODS overlay between successive calls to the sub- 
routine. In the complementary perturbation subroutine indicators that determine 
the logical flow of the program and the time-element array had to be saved before 
each exit and restored at each entry for the same reason. Since the Fortran 
"COMMON" statements could not be used in DODS all subroutine arguments had 
to be passed through the argument list. 

A detailed description of the Brouwer-Lyddane Orbit generator routine 
and the subroutines called by BROWRO follows. The subroutines referenced but 
not described in this document are part of the DODS system and their description 
can be found in the appropriate DODS documentation. 
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DEFINITIVE ORBIT DETERMINATION SYSTEM 


BROWRO — Brouwer- Lyddane Orbit Generation Subroutine 

I. LANGUAGE: 

Fortran IV Level G and Level H 


n. PURPOSE: 

The objective of this subroutine is to compute the position vector (x, y, z) 
and the velocity vector (x, y, z) of an artificial satellite at any time T, when given 
a set of elements at a time T Q called epoch. 

% 

a Q - Semi- major axis 
e 0 - Eccentricity 
- Inclination 
-£ 0 - Mean anomaly 
g 0 — Argument of perigee 
h Q — Longitude of Ascending node 

in. INTERFACE INFORMATION 

A. The Calling Module is ORBGNO 

B. The caUed modules are: 

(1) DRAG — Computes A & drag 

(2) REDUCE — Reduces angle between 0 & 277 

(3) KEPLRI — solves Kepler’s equation 

(4) DATANO — Computes the Double Precision Arctan (y/x) 

(5) PERTFO - Reads complementary Perturbations tape for Brouwer 

(6) MPIDOO — Output intermediate data is requested. 

(7) MPERRO — Error handling routine. 

C. Calling Sequence 

SUBROUTINE BROWRO (EPHEM, K, ELEP, SATID, DPT, ZONALS, 

CF, PLN, IDOBE, IDVICE, IERR, SAVE, INDA). 
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Table I 

Calling Sequence Arguments 
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ELEP(l) 

ELEP(2) 

ELEP(3) 

ELEP(4) 

ELEP(5) 

ELEP(6) 

ELEP(7) 


T 0 Epoch Time in DUT 
a 0 In DUL 

e o 

i0 1 

® 0 v In Radians 

h o f 


DRAG PARAMETERS ARRAY 


DPT(l) - DPT (20) } t 0 , tj ..... t 1# DUT 

DPT(2 1) - DPT (40) } N 2 0 , N 2 j N 2 lg Rad/DUT 2 

DPT (41) - DPT (60) } N ’ 0 , N^ ..... Nj jg Rad/DUT 3 

ZONALS(l) = C 20 2nd zonal harmonic from Harmonic Coeff. file 
ZONALS(2) = C 3 0 3rd zonal harmonic from Harmonic Coeff. file 
ZONALS(3) = C 4 0 4th zonal harmonic from Harmonic Coeff. file 
ZONALS(4) = C 5 0 5th zonal harmonic from Harmonic Coeff. file 


PLN — Perturbation Tape Logical Number 


> 0 Read Pert Tape From Unit Number PLN 
< 0 Do not read Pert Tape 
= 0 Error 


IDOBE(N) = Intermediate data output (IDO) indicator 
IDOBE(N) = 0,No IDO 

IDOBE(N) = i. Output data at every ith iteration 


N = 1 Brouwer elements with secular, long period and short period terms, 
e, i, g, h, 1 

N = 2 Brouwer elements with secular terms only. 

a", e”, 1", g", h", 1" 

N = 3 (L* + G* + H') 

N = 4 L, G, H 

N = 5 Contributions from secular, long, and short period terms 

Se, Si, (Sin i"/2) Sh, e" U 

N = 6 Complementary Perturbations at Request time 

Aa, Ae, Ai, Ag, Ah, A^t 
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N = 7 True Anomaly (f); Eccentric Anomaly (E) at Request time 
N = 8 K2, K3, K4, K5 

NOTE: (1) All IDO will have time as an output. 

(2) IDO's will be Long Format 

(3) [ i] will be constant for all outputs, i.e., [i] will not vary if more 
than one IDO requested. 

IDVICE(l) - Sysout logical file number 
> o use sysout 
< 0 do not use 
= 0 error 

IDVICE(2) - Dedicated printer logical file number. (Treat value same as 
above) 

IDVICE(3) - Tape logical file number (treat value same as above). 
SAVE(l) - SAVE (20) Save area 
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Table 2 

Logical Record Format Working Constants Pool 
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Table 2— Continued 
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INTERFACE BLOCK DIAGRAM 
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FUNCTIONAL ANALYSIS 
BROUWER- LYDDANE FLOW CHART 



COMPUTE: 

’^5 

K 2t K 3 ,K4,K5 


COMPUTE: 

FRACTIONS 


COUNTER FOR 
IDO'S 
= 0 


K-MULTIPLIER 
FOR A-/ orag 


m 
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m 













SEARCH PERT 
TAPE FOR 
ELEMENTS AT 

tn 


compute: 
Ao, Ae, Ai 
A i, Ag, Ah 


COMPUTE 



COMPUTE 
MEAN MOTION 
N USING a" 



COMPUTE 

r z ,r z ,r A j\ 


COMPUTE 

e, e z ,e A 
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8je=(^+g+h} s 0 

ebJ=h i = siN-£ 8h=o 


ID8 = 164 
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FORMULATION 
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I. ZONAL HARMONICS DETERMINATION 

J 2 = “^ 2,0 v/ ^5~; K 2 = 1/2 J 2 Rg 

J3 = “^ 3,0 v ^" : R 3 = “ J3 R e 
J 4 = -C 4>0 /9: K 4 = - 3/ 8 J 4 

Js = -c 5 ,o K 5 = “ Js R l 


where C. 0 is the normalized form of Harmonic Coefficients taken from the 
DODS Harmonic Coefficients Array. 

R e = Radius of earth taken from Constants File 

H. CALL DRAG TO COMPUTE AT TIME T. 

JJKAG 




DRAG “ 


m 3 

T. Z N -> (t - v p 

q ■ 0 ps 2 


where m = 0, 1, 2, • • • 19 
m. COMPUTE MEAN MOTION 

1. Without pert tape N 0 = / ^ 

r a o 

2. With pert tape N = IjL- where a! = a n + Aa/2 

V«) 3 

IV. COMPUTE: 

77 = /l - e* 2 Y 4 = K 4 /a 4 

7 2 =K 2 /a' 2 y; = y 4 /r? 8 

y 2 = y 2 /v 4 y 5 = K s /a* s 
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7 3 =V a ' 3 

y 3 = v ' 7 ? 6 y"s = y s /r > 10 

6 - Cos i" 


|(30 2 - S) + y 2 [f? 2 (- 96 7, + 30 - 90 r, 2 ) 

+ (16 7, + 25t f - 15) + 6 A (144 77 + 25 t? 2 + 105) ]J 
+ e’ 2 y* (jpj (3 + 35 Q A - 30 0 2 )j 
g = N jiL y 4 [[0 2 (126 V 2 - 27) + 6* (385 - 189 t? 2 )] 

- 9 77 2 + 2l] + y\ ^ y 2 [ 9 4 (45 T7 2 + 360 77 + 385) 

+ 6 2 (90 - 192 77 - 126 77 2 ) + (24 77 + 25 tt 2 - 35)] 

+ §(5* 2 -l)]} 

h = N (I) 0(3-75 2 ) (5 - 3 77 2 ) + r' 



^ = 77 N < yl 


x [0 (12 77 + 9 T7 2 - 5) - 9 3 (5 T7 2 + 36 77 + 35)] 
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V. Critical inclination test quantity (I c ) is compared to K c , where K c = .01 


I c = 25 e s y 2 e* 2 / (1 - 5 d 2 ) 2 . 


If I > K then 

C C 


Sj e = {V + g' + h') = e" hi = 8j i = 0 


( Sin f) 


Sh = o 


If I c 1 K c then compute 


S, e, V + g' +h', e" Si, 


>i i, (sini-^ 


8 h 


VI. COMPUTE A1-A10 


A 1 = (|r 2 ' v 2 ^ {1 - 11 e 2 - [(40 0 4 )/(i - 5 d 2 )] } 
/ c\ ( y* \ 

a 2 = hpH — ) v 2 {i - [8 sV(i - 5 e 2 )] - 3 e 2 } 

\12/ \ v / 

\ X 2 / 


7 , 


A 3 = _ (3 e* 2 + 4) 


A, = ^— (3 e" + 4)( {l. 24 S ‘ .»*! 

(i-se 2 ) 


A- (24g4) 

y' \ (1-5 e 2 ) J 
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A io = Sin (i") 1 7 2 

A 6 = (y's/vi) a/4) 


A 


8 “ 



16 £ . 5 * 

(1 - 5 ^) 


COMPUTE B13-B15 


B 13 = e" (A t - A 2 ) 
B 14 =A 7+ (5/64 )(A s )A 10 
Bj 5 = Ag (A 10 ) (35/384) 


COMPUTE A11-A27 


A n = 2 e ' 2 
A 12 = 3 e" 2 + 2 

\ 3 = &2 ( A n) 

A 14 = (5 e' 2 + 2) [0 4 /( 1 - 5 0 2 )] 

A 17 = 8 \/{ 1 - 5 0 2 ) 2 
A 1S = <9 6 e* 2 / (1 - 5 6> 2 ) 2 
A 16 = S 2 /(l - 5 9 2 ) 

A 18 = e " Sin (!') 

A ig = A x8 // (1 * 9 ) 
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A 20 = (1 + 8) Sin i* 

A 21 = e" e 

a 22 = e' 2 e 
A 23 = (A 21 ) Tan 
A 24 - e“ 7) 2 Sin i* 

A 25 = A 12 + 2 

A 26 = 16 ( A ie) + 40 ( A i 7 ) + 3 

A 27 = A 22 ( 1/8 ) Cl1 + 200 ( A l?) + 80 ( A 16 )] 


COMPUTE B1-B12 


B i = V (A, - A 2 ) - |[A n - 400 (A 1S ) - 40 (A 14 ) 

- 11 < A i3>^ (jg) til + 200 (A 1? ) + 80 (A 16 )] (A 22 ) (f 
x y 2 + j[- 80 (A 1S ) - 8 (A 14 ) - 3 (A 13 ) + A n ] (i.) 

+ 12* t A 22^l’ { 

J V ^2 / 

B 2 = t A g) ( A ig) (2 + Tj - e 2 ) + ^ (Ag) (A 19 ) tj 2 

- (H) A 4 ( A is) V 3 + (^) A s + A 6 (A 21 ) Tan (i- 

+ (9 e' 2 + 26) (jjj A 4 (A 18 ) + H (A 3 ) A 21 (A 26 ) 

x Sin i* (1 - 0) 
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B 3 = j[80 (A 1? ) +5 + 32 (A 16 )] (A 22 ) Sin i' ( 6 - 1) 
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B 
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[80 (A 17 ) + 32 (A 16 ) + 5] 


(A 22 ) e* Sin 2 i" 


X 




^8 [^21 ) 


35 

1152 


vn. COMPUTE DOUBLE PRIMED ELEMENTS 


^ =Ut- 1 0 > + n 0 (t - 1 0 ) + + A-e DRAG 

g' = g (t - 1 0 > + g 0 
h* = h (t - t 0 ) + h 0 


f* 


Arctan 


Sin r 
. Cos f* 


where 


Sin f" = 7 ) Sin E* 
Cos f* = Cos E* - e 


r* = a* (1 - e’ Cos E*) 

E" is double primed eccentric anomaly computed from Kepler's equation. 
If Pert tape is being used, 


a 

e 

i 


M 


N 


'0 + A p a 


e o + A p e 
l o + A P 1 



If Pert tape is not used, 


a 


o 


e 


o 


^0 + A p l 
8 0 + A P g 
h 0 + A p h 
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COMPUTE A (SEMI-MAJOR AXIS): 



e = B14 Sin g" + B13 Cos 2g" - B15Sin3g" 

COMPUTE V + g' + h': 

V + g' + h' = l" + g" + h" + B3 Cos 3g" + B1 Sin 2g" + B2 Cos g" 
e" S'C = B4 Sin 2g“ - B5Cosg" + B6 Cos 3g" 



S! = SSini" je" Cos (3f " + 2g") + 3 [e" Cos (2g" + f") + Cos (2g" + 2f ")]} 


(B8 Sin g" + B7 Cos 2g" - B9 Sin 3g" ) 


I") Sh = pr]'<4 B12Cos 3g" + Bll Cos g" + BIO Sin 2g" - \ 

V 1 Cost l L L 


y ' 6 Sin i " 


' / 

X |e(e" sin f -V + f " ) - [3(Sin (2g" + 2f") +e" Sin(2g" + f")) + e" Sin ( 3 f " + 2g" 
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COMPUTE l 

% + g + h = 

x ‘ 

+ 

X 

Se = Sj e + 
X {3e" 



+ [ze 2 
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+ g + h 


r 


^ (77 + 1) 4 e *2 772 

3(1 -e 2 ) s 

Vw 

_ - 


Sin (3f " + 2g" ) 


."V , 

<3 + l r "/ t? 2 + + Sin (2g" + f 


■>(-(.■( 3*4 


2 Sin f " (; 30 2 -l) (v 2 (-p)’ + (-p) + 1^ 


•+7,' 


( 4 ) [{-20 -1 -56*) 


(e* Sin f + f" -V)] + .(3 + 20- 56> 2 ) jy 2 ' ^ [ e - Sin (3f" + 2g") 


+ 3 [sin (2g" + 2f " ) +e" Sin(2g" + f*)l 


r 


4 V 2 


V. I— 




Cos ( 2g* + 2f") 


Cos 2 f" +3 Cos f" + e" 2 Cos 3 f" +e"} 


(l - 6 2 ) [3 Cos (2g"+ f ") + Cos (3f " + 2g")]j. 




+ 3e* Cos 2 f" + 3 Cos f* + e" 2 Cos 3 f 


J 
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COMPUTE e (ECCENTRICITY) 


r (e" U ) 2 + (e* + Se) 


COMPUTE i (INCLINATION) 



If e ^ 0, compute -£ (mean anomaly) 


l 


Arctan 


'e' Si Cos l" + (e" + Se)Sin lA 
(e" + Se ) Cos l" -e" Si Sin l") 


If i = 0 set h = 0 

If i ^ 0 compute h (longitude of the ascending node) 



COMPUTE g (ARGUMENT OF PERIGEE) 

g “ (l + g+h)-l-h 

COMPUTE E (ECCENTRIC ANOMALY) USING KEPLER’S EQUATION. 
COMPUTE f (TRUE ANOMALY) 

f - Arctan 

COMPUTE r (RADIUS VECTOR) 

r 

-g ~ 1 - e Cos E 
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COMPUTE POSITION VECTOR (x, y, z) 

x = r [Cos h Cos (g + f) - Sin h Cos i Sin (g + f)] 

y = r [Cos h Cos i Sin (g + f) + Sin h Cos (g + f ) ] 
z = r [Sin i Sin (g + f)] 

COMPUTE VELOCITY VECTOR (x , y, z) 

x = Sin E (e /Til/ r) [Cos h Cos (g + f) - Sin h Cos i 

x Sin (g + f)] - /l - e 2 /a/x/r [Sin h Cos i Cos (g + f) 

+ Cos h Sin (g + f )] 

y = e Sin E (/a/Vr) [Cos h Cos i Sin (g + f) 

+ Sin h Cos (g + f)] - A - e 2 (/ajl/r) [Sin h Sin (g + f) 
- Cos h Cos i Cos (g + f)] 
z = e Sin E (/ a/x/r) Sin i Sin (g + f) 

+ [A - e 2 (/ aJI/r ) Sin i Cos (g + f)] 
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RESTRICTIONS AND LIMITATIONS 

I. Mathematical Restrictions: None 

II. Data Restrictions: None 

III. Hardware Restrictions: None 

IV. Programming Language Restrictions: None 
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DEFINITIVE ORBIT DETERMINATION SYSTEM 
PERTFO - Complementary Perturbations Tape Read Routine 
I. LANGUAGE: 

Fortran IV, Level G and Level H 
n. PURPOSE: 

This subroutine reads the complementary perturbations tape for the 
Brouwer- Lyddane Orbit generator. 

HI. INTERFACE INFORMATION 

A. Calling Module is BROWRO 

B. Called Modules are 

(1) INTPLO — performs backward difference interpolation for Comple- 
mentary Perturbations Tape Read Routine 

(2) MPERRO — Error handling routine from DODS 

(3) TCONVO — Time Conversion Routine from DODS 

C. Calling sequence 

Subroutine PERTFO (PLN, SATID, TIME, KMULT, B, IERR, EPHEM) 
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DEFINITION OF ARRAYS: 
B ARRAY 


B (1) = a in DUL 

P 


B (2) = e p 
B (3) = i p (Rad) 
B (4) = l (Rad) 
B (5) = g D (Rad) 


> Elements from perturbations tape. 


B (6) = h p (Rad) 


PLN — Perturbations tape unit number 

PLN > 0 Read Pert tape on this unit 
PLN < 0 Do not read Pert tape 
PLN = 0 Error 
EPHEM ARRAY 


Array to store variables that might be destroyed by DODS overlay. 
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PERTFO FLOWCHART 



m 
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POSITION 
TAPE FOR 
FIRST TIME 
SEARCH 


SET 

INDICATOR 
190 =0 



RETURN 





f REQ > 
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READ 6 RECORDS 
FROM PERT TAPE 
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BROUWER PERTAPE FORMAT 


Word 

No. 

Contents 

i 

Record #1 - Header Record 

1 

Fortran Word Count 

2 

Time Increment - days 

3 

Month 

4 

Day 

5 

Year 

6 

Satellite ID number 

7 

Input Semi-major Axis - E.R. 

8 

Input Eccentricity 

9 

Input Inclination - degrees 

10 

Input right ascension of the ascending node - degrees 

11 

Input argument of perigee - degrees 

12 

Input mean anomaly - degrees 

13 

Input time from midnight - days 

14 

Input period - minutes 

15 

No. of records on tape excluding header and trailer 

16 

Delta mean anomaly option indicator 


(1 - delta drag mean anomaly not computed on tape 


0 - delta drag mean anomaly computed on tape.) 

Record 2 to N 

1 

Fortran Word count 

2 

Time in seconds from epoch 

3 

A (semi-major axis) - ER 

4 

e (eccentricity) 

5 

i (Inclination) - 77/2 - i = w/2 

6 

AM (mean anomaly change from t 0 ) 0 = AM = 277 

7 

co (argument of perigee) 0 = co= 2-n 

8 

G (right ascension of the ascending node) 0 ^ fi ^ 2n 


Last Record - Trailer Record 

1 

Fortran Word count 

2 

.99999999 x 10 30 

3-8 

Dummy Words 


39 













PAGE: 12-12 
DATE: 12 January, 1970 
SYMBOL: PERTFO 
CONTRIBUTOR: E. A. Galbreath 
GSFC 


RESTRICTIONS AND LIMITATIONS 

l. Mathematical Restrictions: 

None 

n. Data Restrictions: 

There must be at least 3 data records on the perturbation tape before the 
start time for the run and at least 3 data records on tape after the end 
time for the run. 

m. Hardware Restrictions: 

None 

IV. Programming Language Restrictions: 

The perturbation tape may be backspaced during the program, therefore the 
pert tape must be unblocked. Blocked tapes give undetermined results if 
backspaced. 
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DEFINITIVE ORBIT DETERMINATION SYSTEM 
INTPLO - Backward Difference Interpolation Function 

I. LANGUAGE: 

Fortran IV, Level G and Level H 

II. PURPOSE: 

This subroutine interpolates for the elements a ,e , i , -t ,g,h when 

r P ' P ’ P w P °P* P 

given a request time (t REQ ) between two times on the tape. 

III. INTERFACE INFORMATION 

A. Calling Module is PERTFO 

B. Calling Sequence SUBROUTINE INTPLO (TIME, A, B, TSUBO, DELTA) 
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DEFINITION OF ARRAYS: A(I, J) 



B ARRAY 


B(l) = a 
B(2) = e 

B(3) = i l 
B(4) = l 
B(5) = g p 
B(6) = h 

p 


in DUL 


Radians 
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FUNCTIONAL ANALYSIS 

Backward Difference Interpolation 



1 3 V times from the perturbation tape 





f (t 5 ) 

f(t 4 ) 

f(ts) 

f(t 2 ) 

f(tx) 


> 




elements from the perturbation tape at time t. . 


t=^ request time 

A t =^> time increment from perturbation tape 
A k f (tj)=*the k th difference of the elements A k_l f (t . + - A k-1 f (t .) 

where A 0 f (t.) s f (t.) 


DIFFERENCE TABLE 


f (t 6 )> 


f(t$) \ ^>A2f(t 4 K 

^>Af(t 4 )^ \A 3 f(t 3 ) 

^A2f(t 3 )^ ^A 4 f(t 2 X 

^>Af(t 3 )<^ J>A 3 f(t 2 )^ 

f(t3 < \A*f(t 2 )<^ >A 4 f(t 

\Af(t 2 )^ J)A 3 f( tl ) 

^Af (t 2 )^ 
f (t x K 


' A 5 f (t t ) 
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(element) p = f (t 6 ) + Af (t 5 ) (tt 0 ) + A 2 f (t 4 ) 


(tt 0 ) (tt t ) 
2 


+ A 3 f (t 3 ) 


(tt 0 ) (tt x ) (tt 2 ) 


+ A 4 f (t 2 ) 


(tt 0 ) (ttj) (tt 2 ) (tt 3 ) 


24 


+ A 5 f (tj) 


(tt 0 ) (ttj) (tt 2 ) (tt 3 ) (tt 4 ) 


120 


Where 


t * 0 = t-t 6 

ttj = At + tt 0 
tt 2 = At + ttj 
tt 3 = At + tt 2 
tt 4 = At +tt 3 


RESTRICTION AND LIMITATIONS 

I. Mathematical Restrictions: None 

II. Data Restrictions: None 

III. Hardware Restrictions: None 

IV. Programming Language Restrictions: None. 
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DEFINITIVE ORBIT DETERMINATION SYSTEM 
Drag - AL Drag Subroutine 


I. Language: 

Fortran IV, Level G and Level H 

II. Purpose: 

This subroutine computes A^ drag for time t = t - t 0 
from input parameters t p q and N p . 

III. Interface Information 

A. Calling Module is BROWRO 

B. Called Module is REDUCE which reduces angle between o & 2n 

C. Calling Sequence 

Subroutine Drag (DPT, PI2, DRAGL, TO, T, KMULT) 
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DEFINITION OF ARRAYS 


DPT (1) 
DPT (20) 


> t 0 , t x , . • • , t ig 


DPT (21) | 

| > N 2 o . N 2 j , . . . , N 2 i 9 


DPT (40) 


DPT (41) 


^ 3 , 0 ’ ^ 3 , 1 ’ ’ ' ' ’ ^3 , 


19 


DPT (60) 
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INTERFACE BLOCK DIAGRAM 
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DRAG FLOWCHART 
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PROCESSED 


RETURN 
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FUNCTIONAL ANALYSIS 

I. Formula for computing 

m 3 

A ^drag = ZZ ZZ N p-q^ t “ t q^ P 

q = 0 p= 2 

where m = 0, 1, 2, • • • • 19 


RESTRICTIONS AND LIMITATIONS 

I. Mathematical Restrictions: None 

II. Data Restrictions: None 

III. Hardware Restrictions: None 

IV. Programming Language Restrictions: None 
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